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1.  Introduction 

This  is  the  final  report  documenting  accomplishments  and  progress  as  a  result  of  research 
performed  under  the  AFOSR  Grant  FA9550- 12- 1-0461,  over  a  period  of  50  months,  spanning 
from  30  September  2012  through  30  November  2016,  following  previous  and  initial  guidance  by 
Dr.  Julian  Tishkoff  and  subsequently  by  Dr.  Chiping  Li  serving  as  program  managers.  The 
research  and  work  represented  a  closely  knit  collaborative  effort  between  the  California  Institute 
of  Technology,  with  Prof.  Paul  E.  Dimotakis  serving  as  Principal  Investigator  (PI),  and  the 
University  of  Minnesota,  with  Prof.  Graham  V.  Candler  serving  as  Co-Principal  Investigator 
(Co-PI). 

The  research  effort  focused  on  high-speed  air-breathing  propulsion  and  turbulent-combustion 
flows  and  environments,  and  was  comprised  of  several  parts  in  which  significant 
accomplishments  were  made: 

1.  An  experimental  effort  focusing  on  investigations  in: 

a.  supersonic-flow  turbulent  mixing,  in  both  chemically  reacting  and  non-reacting 
flow,  as  well  as 

b.  small-scale  fully  resolved  scalar  mixing  experiments  aimed  at  providing  guidance 
to  and  validation  of  subgrid-scale  model  of  scalar  transport  and  mixing. 

2.  An  equipment  and  instrumentation  effort  comprised  of  experimental  technology  and 
high-speed  data-acquisition  developments,  as  well  as  modeling  and  simulation  advances 
that  permitted  new  validation  approaches  of  large-eddy  simulation  (LES)  numerical 
modeling  based  on  detailed  forward-modeling  of  shadowgraph-schlieren  optical  imaging. 

3.  A  theoretical  and  numerical-simulation  effort  focused  on  subgrid-scale  (SGS)  model  and 
LES  developments: 

a.  non-reacting  variable-density  compressible  and  incompressible  flows; 

b.  subgrid-scale  modeling  of  turbulent  combustion,  extending  to  hydrocarbon 
combustion,  and  yielding  a  new  framework  for  such  modeling  dubbed  the 
evolution-variable  manifold  (EVM)  method; 

c.  development  of  numerical  flux  functions  and  implementation  of  the  evolution- 
variable  manifold  (EVM)  approach  in  the  University  of  Minnesota  US3D 
computational  fluid  dynamics  code; 

d.  Simulations  of  low-speed  and  supersonic  Caltech  hydrogen-fluorine  reacting 
mixing-layer  experiments 

e.  Development  of  a  novel  SGS  modeling  approach  for  LES  modeling  of 
compressible  and  variable-density  flows 

4.  Theoretical  and  numerical  investigations  of  unphysical  scalar  excursions  that  arise  as  the 
result  of  dispersive  errors  encountered  in  the  transport  of  passive  and  active  scalars  in 
finite-difference  numerical-simulation  schemes,  and  preliminary  investigations  of 
mitigation  measures  for  both  passive  and  active  scalars. 
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5.  Large-scale  direct  numerical  simulation  (DNS)  modeling  of  turbulent  flows  that  mix  two 
fluids  spanning  density  ratios,  R  =  p1/p2,  in  the  range  of  1.005  <  R  <  10,  based  on 
spectral  numerical  methods  that  are  free  of  the  dispersive  errors  outlined  under  Item  4, 
above.  Turbulent  hydrocarbon  combustion,  for  example,  entails  flow  spanning  density 
variations  of  order  R  =  punburnt/ Pburnt  —  7  for  reactions  between  species  starting  at 
normal  temperatures  and  pressures.  These  investigations  aim  to  probe  and  understand 
dynamics  such  as  baroclinic  torques  not  encountered  in  (near-)uniform-density  flows. 

To  date,  this  work  has  been  documented  in  papers  in  various  stages,  i.e.,  published,  submitted,  or 
in  preparation  at  this  writing,  as  well  as  in  topical  conference  presentations  and  papers. 

Research  completed  under  this  project  enables  several  advances.  Importantly,  efficient  LES 
modeling  of  complex-geometry  high-speed-propulsion  flow  paths  with  turbulent  combustion  in 
scramjet-type  combustors  is  now  possible.  The  EVM-LES  approach  could  in  the  future  be 
extended  to  long-chain  hydrocarbon  fuels  using  recent  advances  in  the  modeling  of  fuel 
breakdown  and  pyrolysis,  followed  by  oxidation  of  short-chain  species.  Crucially,  the  present 
research  developed  an  accurate  approach  for  correctly  linking  multi-dimensional  tabulated  EVM 
data  to  thermodynamic  processes  that  must  be  correctly  captured  by  and  represented  in  the 
numerical  method. 

The  effort  benefitted  from  progress  under  previous  support  by  AFOSR  Grant  FA9550-07-1-0091 
and  the  AFOSR  DURIP  Grant  FA9550-10-1-0553,  as  well  as  use  of  equipment  and 
instrumentation  previously  developed  under  a  NSF  Major  Research  Instrumentation  Grant  EIA- 
0079871.  More  recently  and  concurrently,  it  benefitted  by  support  under  the  DOE  Research 
Grant  Award  DE-NA0002382  that  continues  at  this  writing  and  support  for  I.  Gat  at  Caltech  by 
the  NSF  Graduate  Research  Fellowship  Program  Grant  Award  DGE- 1144469.  Further  support 
was  provided  by  the  Caltech  academic  program,  and  the  Caltech  Northrop  Chair  in  Aeronautics, 
as  well  as  through  the  McKnight  Presidential  Professorship  at  the  University  of  Minnesota. 

Computing  support  was  also  provided  by  the  Blue  Waters  sustained-petascale  computing  project, 
supported  by  NSF  Awards  00-0725070  and  ACI- 1238993,  and  the  state  of  Illinois.  Blue 
Waters  is  a  joint  effort  of  the  University  of  Illinois  at  Urbana-Champaign  and  its  National  Center 
for  Supercomputing  Applications.  Simulations  performed  on  Blue  Waters  were  under  NSF 
PRAC  award  number  ACI-1440083.  Computations  were  also  performed  on  the  Caltech  Zwicky 
computer  cluster,  supported  by  NSF  MRI-R2  Award  PHY-0060291  and  by  the  Sherman 
Fairchild  Foundation.  The  work  was  also  supported  by  the  Cray  Trinity  system  of  the  Alliance 
for  Computing  at  Extreme  Scale  (ACES),  a  partnership  between  Los  Alamos  National 
Laboratory  and  Sandia  National  Laboratories  for  the  U.S.  Dept,  of  Energy's  NNSA.  Data  storage 
and  visualization  were  facilitated  by  a  computer  cluster  integrated  by  Daniel  Lang  and  developed 
through  support  by  NSF  MRI  Grant  EIA-0079871,  AFOSR  DURIP  Grant  FA9550-10-1-0553, 
and  support  by  the  AFOSR  and  DOE  grants  mentioned  above.  Computational  support  was  also 
provided  by  the  University  of  Minnesota  Supercomputing  Institute,  the  University  of  Minnesota 
Aerospace  Engineering  &  Mechanics  computer  clusters,  and  the  DoD  HPC  systems  through  the 
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AFOSR  grant  support.  We  would  finally  like  to  acknowledge  discussions  with  Profs.  Dan 
Meiron  and  Dale  Pullin,  and  a  collaboration  in  the  computations  with  Prof.  Christian  Ott. 

The  research  teams  comprised  of  faculty,  students,  post-docs,  and  other  researchers  and 
collaborators  are  grateful  to  the  AFOSR  and  other  support  provided  that  enabled  the  important 
progress  we  report  on  below. 
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2.  Research  and  progress  under  this  grant 

Research,  accomplishments,  and  progress  realized  under  this  Grant,  as  outlined  in  the 
Introduction,  are  documented  below. 

2.1  Experimental  effort 

This  section  summarizes  the  experimental  effort  performed  under  this  Grant  and  documents  its 
results. 

2.1.1  Experiments  in  supersonic  flow  and  mixing 

Several  experiments  were  undertaken  in  supersonic  flow  investigating  the  behavior  and  of  both 
reacting  and  non-reacting  flows.  These  were  primarily  focused  on  transverse-jet  injection  in 
supersonic  cross-flows;  a  rich  flow  environment  relevant  to  scramjet  propulsion  that 
encompasses  many  of  the  dynamics  investigated  separately,  either  experimentally  or 
computationally,  as  discussed  below.  See  Figure  1. 


Figure  1.  Transverse-jet  injection  in  supersonic  crossflow  with  an  inflow  Mach  number  of  M}  —  1.5. 


These  experiments  necessitated  modifications  to  the  test-section  to  enable  a  modular 
arrangement  for  interchangeable  blocks  for  jet  injection  at  90°,  or  inclined  at  30°,  with  round  jets 
of  different  diameters,  specially  designed  jet  inlets  to  avoid  complications  from  sharp  inlet  edges, 
diamond  jet  arrays,  or  single  jets,  or  jets  in  combination  (Cymbalist  2016).  They  also  required 
the  development  of  new  experimental  technologies  and  instrumentation  that  were  undertaken 
under  support  of  this  grant  as  well  as  support  from  the  AFOSR  DURIP  Grant  FA9550-10-1-0553 
and  developments  under  previous  AFOSR  and  NSF  support,  as  outlined  in  the  Introduction. 

Optical  diagnostics  employed  combined  shadowgraph-schlieren  techniques  with  high-speed/- 
resolution  digital-image  recording,  with  concurrent  chemiluminescence  imaging  recorded  on  a 
gated  image-intensified  camera.  These  diagnostics  exploited  counter-propagating  optical  beams 
with  a  west-to-east  propagation  for  the  shadowgraph-schlieren  data  and  a  west-to-east 
propagation  for  the  chemiluminescent  data,  further  separated  from  each  other  exploiting  optical 
wavelength  differences  using  dichroic  optics.  A  schematic  diagram  of  the  combined  optical 
system  is  depicted  in  Figure  2  from  Cymbalist  (2016). 
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Figure  2.  Schematic  of  combined  shadowgraph-schlieren  and  chemiluminescence  optical  diagnostic 
system  (Cymbalist  2016,  Fig.  2.4). 

Images  were  co-registered  on  the  computer  exploiting  fiducials  in  the  images  and  could  either  be 
superimposed  or  processed  separately  to  yield  instantaneous  (/rs)  density  fields,  or 
chemiluminescence  fields,  as  illustrated  in  the  transverse-jet  injection  in  a  supersonic  (M,  =  1.5) 
crossflow  in  Figure  3.  The  Mach-wave,  bow-shock,  acoustic-wave  system,  and 
chemiluminescence  intensity  distribution  are  discernible  in  that  figure. 

Three  distinct  regions  are  discernible  in  the  bottom  annotated  chemiluminescence  frame- 
averaged  image,  with  { m  a  spatial-delay  length.  Also  interesting  in  Figure  3  is  the  near-wall 
reaction  zone  that  can  be  seen  to  even  extend  slightly  upstream  of  the  jet,  marking  the  horseshoe 
recirculation  zone  around  the  jet  and  the  long  residence  times  giving  reactions  the  opportunity  to 
take  place  in  that  region. 

Chemical-kinetics  as  well  as  mixing  scaling  analyses  were  explored  as  causes  to  investigate  the 
dependence  of  the  delay  length  as  this  is  an  important  parameter  in  scramjet  combustion  and 
propulsion. 

The  hydrogen-fluorine  chemical  system  is  relatively  simple  and  kinetically  fast,  even  at  low 
reactant  concentrations,  so  it  was  possible  to  exclude  kinetics  by  performing  kinetics-delay 
estimates  as  the  cause  for  this  spatial  delay.  This  was  further  confirmed  by  changing  kinetics 
rates  experimentally  by  increasing  reactant  concentrations,  which  had  a  very  small  effect. 
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Figure  3.  Transverse-jet  injection  into  a  supersonic  N2  crossflow  stream  with  an  inflow  Mach  number  of 
Mi  =  1.5.  Wj  denotes  the  jet-fluid  (He)  molar  mass  while  dj  —  0.20"  is  the  jet  diameter.  Top:  Single 
shadowgraph-schlieren  digital  image  with  a  superimposed  chemiluminescence  image  (from  Cymbalist 
2016,  Fig.  2.15).  Bottom:  Frame-averaged  chemiluminescence  image  from  the  same  run  (from  Cymbalist 
2016,  Fig.  3.1).  Color  legend  (bottom,  right)  registers  intensity. 


A  separate  investigation  explored  the  possibility  that  the  spatial  delay  is  attributable  to  a  mixing 
delay,  which  the  chemical  reaction  must  await  before  any  significant  reaction  occurs.  This  was 
confirmed  by  an  experimental  investigation  that  varied  the  jet  Reynolds  number, 


Rej 


P j  P]  dj 
Pi 


(2.1.1) 


where  pj  is  the  jet-fluid  density,  Uj  is  the  jet  velocity,  dj  is  the  jet  diameter,  and  pj  is  the  jet 
(shear)  viscosity,  by  varying  the  jet  diameter  as  well  as  the  jet-fluid  density  using  a  different 
mean  molar  mass  for  the  jet  fluid. 


Figure  4  plots  the  delay  length  scaled  by  the  (geometrical)  jet  diameter,  Tm/dj,  vs.  a  range  in  jet 
Reynolds  number,  in  terms  of  log(Rej).  As  can  be  seen,  while  the  general  trend  suggests  a 
Reynolds-number  dependence  and  a  fluid-mechanics  cause  for  this  delay,  the  plot  does  not 
indicate  good  collapse  with  respect  to  these  parameters. 


If  the  spatial  delay  is  the  result  of  turbulent  mixing  that  occurs  in  the  far  field,  however,  the 
proper  scaling  length  that  also  takes  into  account  that  the  mean  density  in  the  jet  becomes 
dominated  by  that  of  the  entrained  fluid  is  the  jet-source  diameter,  i.e.. 


2rhj 

y/ltPaoJ] 


(2.1.2) 
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where  rhj  is  the  jet-fluid  mass  flux  through  the  jet  nozzle,  pm  is  the  density  of  the  far-field  fluid 
that  the  jet  entrains,  and  /j  is  the  jet  momentum  flux  (thrust).  These  are  integral  quantities  felt  in 
the  jet  far  field  and  parameterize  entrainment  and  mixing  far  from  the  jet  nozzle.  This  scaling  is 
found  to  correlate  jet  flame  lengths  ranging  from  gas -phase  reactions,  to  hydrocarbon 
combustion,  to  reacting  jets  in  water.  See  Dimotakis  (2000)  and  references  therein. 

Figure  5  plots  the  ratio  of  the  delay  length  to  d*,  vs.  log(f?ej),  indicating  very  good  scaling  in 
those  terms  and  basically  confirming  that  the  spatial  delay  is  attributable  to  turbulent  mixing,  as 
chemical  reactions  must  await  that  process  before  proceeding. 


Figure  4.  Mixing-delay  length  normalized  in  terms  of  the  physical  jet-orifice  diameter,  dj,  vs.  the  jet 
Reynolds  number  (Eq.  1).  Bars  reflect  uncertainty  based  on  the  variation  of  £m. 

Jet-inclination  effects  were  studied  by  also  investigating  and  comparing  the  behavior  of  normal 
jet  injection  (0j  =  90°)  with  that  of  jets  inclined  at  an  angle  of  0j  =  30°.  Figure  6  compares  the 
two  cases  for  reacting  jets,  superimposing  shadowgraph-schlieren  data  with  chemiluminescence 
data  and  demonstrates  the  large  difference  in  the  cross-flow  disturbance  between  normal 
injection  and  that  from  the  inclined  jet. 
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Figure  5.  Variation  of  ^m/d*  vs.  jet  Reynolds  number  (same  data  as  in  Figure  4). 
This  should  help  in  scaling  and  design  studies  of  scramjet  combustor  systems. 
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Figure  6.  Superimposed  shadowgraph-schlieren  and  chemiluminescence  images  of  jet  injection  into  a 
supersonic  stream.  Top:  normal  injection.  Bottom:  30°  inclined  jet.  The  A-shock  structure  (*  symbol)  and 
acoustic  waves  (A  symbols)  are  evident  in  the  normal-jet  case,  but  absent  in  the  inclined-jet  flow. 
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Figure  7  compares  the  corresponding  (enhanced)  chemiluminescence  images  for  the  two  cases 
demonstrating  the  longer  reaction-delay  length  for  the  inclined-jet  case  but  also  the  lower 
penetration  into  the  freestream,  an  important  consideration  for  air-breathing  scramjet  propulsion. 
On  the  other  hand,  the  deeper  penetration  of  normal-jet  injection  is  accompanied  with 
considerable  total-pressure  losses,  as  can  be  inferred  from  the  shock  structure  in  Figure  6. 
Finally,  we  note  that  absent  in  the  inclined-jet  injection  is  a  reaction  zone  upstream  of  the  jet- 
injection  station,  indicating  the  difference  in  near-wall  mixing  and  combustion  between  the  two 
flows. 
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Figure  7.  Chemiluminescence  images  from  chemically  reacting  jets  injected  into  a  Mj  =  1.5  cross-flow. 
Top:  normal  injection  if)]  —  90°,  same  run  as  in  Figure  3).  Bottom:  inclined  injection  at  0j  =  30°. 


2.1.2  Experiments  in  small-scale  mixing 

Experiments  in  small-scale  mixing  continued  under  support  of  this  grant,  albeit  at  low  priority  in 
comparison  to  other  focused  efforts  described  in  this  report.  Results  from  these  experiments  are 
shedding  light  into  the  small-scale  scalar-mixing  structure  providing  guidance  to  subgrid-scale 
(SGS)  model  development. 

The  experiments  are  conducted  in  water  and  rely  on  a  facility  (Figure  8)  developed  for  the 
purpose  to  study  the  simplest  turbulent  flow,  namely,  grid-generated  turbulence,  with  the 
capacity  to  generate  a  large  range  of  Reynolds  numbers  based  on  the  grid-mesh  size,  £g,  i.e., 
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500  <  ReM  = 


Mi 

V 


<  56,000, 


(2.1.3) 


where  Ug  is  the  grid  velocity,  equivalent  to  the  flow  velocity  in  a  water-tunnel  frame,  and  v  is  the 
kinematic  viscosity.  This  Reynolds  number  range  spans  flows  from  below  to  above  the  mixing 
transition  (Dimotakis  2000). 


The  experiments  yield  fully  resolved  3D  and  4D  scalar-field  data  at  the  lower  Reynolds  numbers 
and  a  window  to  processes  typically  beyond  experimental  reach.  Measurements  rely  on  laser- 
induced  fluorescence  of  a  scalar  tracer  injected  just  downstream  of  the  grid  that  then  disperses 
and  mixes  in  the  ensuing  grid-turbulence  field.  Custom-designed  digital-imaging  equipment  and 
techniques  developed  starting  with  support  by  the  NSF  MRI  Grant  EIA-0079871  and  the 
AFOSR  DURIP  award  provided  the  instrumentation.  Full-resolution  measurements  impose 
challenging  data-acquisition  and  high-speed  storage  requirements  derived  from  the  sustained 
measurement  rate  of  2  x  108  measurements  (pixels)  per  second.  While  such  data  rates  are  not 
uncommon  today,  what  remains  uncommon  is  the  need  to  sustain  them  for  many  minutes  to  yield 
the  data  for  the  reconstruction  of  the  required  3D  and  4D  fields.  The  experiments  themselves 
were  conducted  (at  low  priority)  under  the  AFOSR  grant  support  documented  in  this  report. 


Figure  8.  Grid-turbulence  scalar-dispersion  experimental  facility. 

The  facility  optics  developed  allow  both  streamwise  and  transverse  cuts  to  be  made  by  rotating  a 
single  Littrow  prism,  while  retaining  optical  alignment.  The  parabolic  mirror  is  an  astronomical- 
quality  reflector  and  produces  a  perfectly  collimated  beam  obviating  complicated  image 
processing  that  was  attempted  in  the  early  phases  of  this  experiment. 

Figure  9  below  reproduces  a  streamwise  cut  of  the  scalar  field  and  attests  to  the  high  signal-to- 
noise  ratios  and  image  quality  attained.  The  field  of  view  in  this  figure  spans  from,  roughly,  16  to 
24  grid-mesh  lengths  (16  <  x/M  <  24)  downstream.  The  mesh  Reynolds  number  is  a  little 
below  the  mixing  transition.  Evident  in  this  image  is  the  degree  of  unmixedness  in  this  flow,  with 
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completely  unmixed  fluid  crossing  the  full  extent  of  the  downstream  scalar  plume,  as  well  as  the 
very  high  local-instantaneous  scalar-concentration  values. 
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Figure  9.  Scalar  field  vertical  streamwise  cut,  16  <  x/M  <  24  downstream  of  the  grid.  Mesh  Reynolds 
number,  ReM  —  4200.  Color  codes  local  scalar  concentration  (Dimotakis  and  Lang,  unpublished  data). 


Figure  10.  Scalar  field  vertical  cross-stream  cut,  x/M  =  20,  downstream  of  the  grid.  Mesh  Reynolds 
number,  ReM  —  4200.  Color  codes  local  scalar  concentration  (Dimotakis  and  Lang,  unpublished  data). 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


13 


Cross-stream  cuts  were  the  most  convenient  to  assemble  to  form  the  3D  images  that  are  under 
investigation  at  this  writing.  An  example  of  such  a  laser- sheet  slice  is  reproduced  in  Figure  10, 
recorded  at  a  distance  x/M  ^  20  downstream,  mid-distance  along  the  field-of-view  that 
produced  the  data  in  Figure  9,  at  the  same  flow  conditions  that  yielded  the  data  in  Figure  9. 

Results  indicate  a  qualitative  transition  in  small-scale  structure  as  the  flow  crosses  the  mixing 
transition,  as  can  also  generally  be  inferred  from  chemically  reacting  experiments  sensitive  to 
small-scale  mixing. 

This  research  is  continuing  as  an  unsupported  effort  at  this  writing,  albeit  slowly,  and  we  hope  to 
be  able  to  document  the  important  results  in  the  near  future. 

2.2  New  instrumentation,  diagnostics,  and  experimental  techniques 

As  part  of  the  diagnostics-development  effort,  successive  pairs  of  schlieren- shadowgraph 
images,  spaced  by  a  short  interframe  time  interval  (A tIF  =  5  —  6/rs)  apart,  were  correlated  by 
extending  previously  employed  methods  to  infer  the  velocity  fields.  Such  correlations  track  the 
convection  of  index-of-refraction  (IoR)  interfaces  in  the  flow,  which  can  be  discerned  in  Figure 
3,  top,  and  in  Figure  6. 


Figure  11.  A-B  image  pair  recorded  spaced  by  an  interfame  time  interval  of  A  tIF  =  6/<s.  Large-scale 
structures  associated  with  the  injected-jet  flow  convect  (are  approximately  displaced  in  the  streamwise 
direction)  with  little  change  in  shape  during  such  intervals. 

Such  correlation  techniques  have  been  proposed  in  the  past  (e.g.,  Townsend  1936)  and  employed 
in  estimating  convection  velocities  of  flow  structures  in  various  flows,  including  supersonic  flow 
(e.g.,  Papamoschou  1989),  as  well  as  in  flow  estimation  for  jets-in-crossflow  (Gruber  et  al.  1997, 
Ben-Yakar  et  al.,  2006).  These  methods  relied  on  manual/subjective  identification  of  flow 
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structures  used  to  track  convection  velocities.  We  decided  to  develop  an  objective  velocity- 
estimation  method  instead,  as  described  below. 

Part  of  the  challenge  is  that  not  all  features  in  such  images  convect  similarly,  such  as  Mach 
waves  that  are  (almost)  stationary,  boundary-layer  structures  that  move  with  a  fraction  of  the 
freestream  velocity,  etc.  The  technique  developed  more-closely  derives  from  that  of  Jonassen  et 
al.  (2006),  adapted  for  the  special  flows  of  interest  here  and  the  combination  of  flow  structures 
that  should  and  should  not  be  taken  into  account. 

A  typical  A-B  image  pair  is  shown  in  Figure  11,  spaced  by  an  interframe  time  of  At1F  =  6/rs. 
Large-scale  flow  structures  remain  well-correlated  across  such  short  time  intervals,  while  quasi¬ 
stationary  structures  registered  as  (almost)  fixed  in  space  and  can  be  relatively  easily  removed 
from  being  considered  for  processing. 
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Figure  12.  Binning  structure  (top)  used  to  produce  average  convective  velocity  estimates  (bottom). 

The  region  downstream  of  the  jet  is  binned  into  rectangular  regions  over  which  the  schlieren- 
image  correlation  velocimetry  (SICV)  results  are  averaged  to  produce  the  convective-velocity 
estimates.  The  bins  and  typical  results  are  shown  in  Figure  12.  The  freestream  velocity  can  be 
reliably  estimated  from  Mach-wave  angles  and  is  also  indicated  in  the  bottom  panel  of  the  figure. 

An  important  conclusion  is  that  velocities  downstream  of  the  supersonic -jet  body  quickly 
approach  that  of  the  freestream  so  turbulent  mixing  that  ensues  between  the  freestream  and  jet 
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fluids  takes  place  in  an  essentially  incompressible  (subsonic)  flow  environment.  Further  details 
and  discussion  can  be  found  in  Cymbalist  (2016)  and  Cymbalist  et  al.  (2017,  in  preparation). 

Both  shadowgraph-schlieren  and  convection-velocity  experimental  data  were  utilized  in  a  novel 
way  that  relied  on  forward-modeling  of  beam  propagation  through  the  optical  system  (Figure  2) 
based  on  ray-tracing  in  the  geometrical-optics  approximation,  as  well  the  estimated  convective 
velocity  fields,  as  illustrated  in  the  bottom  panel  of  Figure  12. 


Figure  13.  Half-volume  results  of  the  jet-injection  LES  modeling.  Graphic  plots  isosurfaces  of  the  Q- 
criterion  colored  by  streamwise  velocity. 


The  supersonic  jet-injection  flow  was  simulated  using  the  US3D  U.  Minnesota  code  to  produce 
instantaneous  3D  flow  fields,  as  depicted  in  Figure  13  that  shows  half- volume  results  of  jet- 
injection  LES  modeling  plotting  isosurfaces  of  the  Q-criterion  colored  by  streamwise  velocity. 

The  jet-injection  region  has  the  highest  velocity  (red)  and  crossflow  is  into  the  page.  The  center 
mid-span  on  the  right  plots  contours  of  |Vp|  showing  shock  structures  and  eddy  interfaces  in  this 
region. 

Geometrical-optics  forward  modeling  was  then  used  to  propagate  the  optical  beam  through  all 
elements  of  the  beam-  and  image-forming  system,  and  though  the  instantaneous  index-of- 
refraction  (IoR)  field  in  the  turbulence  to  compute  simulated  shadowgraph-schlieren  images. 
Typical  results  by  way  of  illustration  are  shown  in  Figure  14. 

In  practice,  an  image  is  almost  never  a  pure-shadowgraph  or  pure-schlieren  image  and  the 
simulations  also  modeled  how  the  knife-edge  was  set  to  match  the  optical  choices  made. 

An  important  application  of  this  technique  was  processing  pairs  of  images  computed  based  on 
the  LES  fields,  as  above,  spaced  by  the  same  interframe  time  interval,  and  subjecting  them  to  the 
same  image-correlation  velocimetry  technique  described  above.  This  produced  estimates  of 
experimentally  measured  vs.  computed  velocity  fields  independently  of  each  other. 
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Figure  14.  Experimental  and  computational  schlieren  and  shadowgraph  images.  Top  and  upper  middle: 
Experimental  and  computational  shadowgraphs.  Lower-middle  and  bottom:  Experimental  and 
computational  schlieren. 

Typical  results  are  depicted  in  Figure  15  that  allows  a  comparison  of  the  two  estimates.  As  can 
be  seen,  except  for  the  region  just  downstream  of  the  jet  nozzle  where  the  lower  computational 
resolution  did  not  allow  the  same  degree  of  detail  to  be  discerned,  the  results  are  in  good 
agreement.  Also  in  fairly  good  agreement  is  the  jet-penetration  extent.  LES-SGS  modeling 
options  were  not  explored  to  investigate  the  relative  merits  of  various  SGS  model  options. 
Further  details  are  documented  in  Luthman  et  al.  (2017,  in  preparation). 
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Figure  15.  Streamwise  velocity  of  convected  refractive -index  interfaces  overlaid  on  an  instantaneous 
shadowgraph  under  the  same  flow  conditions.  Top:  experiment.  Bottom:computation. 

2.3  Large-scale  direct  numerical  and  large-eddy  simulations 

Part  of  the  goal  of  the  present  research  effort  was  not  only  to  probe  the  behavior  of  flows  with 
attributes  relevant  to  high-speed  air-breathing  propulsion  applications,  but  also  to  investigate 
generic  fundamental  behavior  of  flows  that  entail  variations  in  density  that  cannot  be  ignored,  as 
well  as  suitable  subgrid-scale  (SGS)  models  for  use  in  large-eddy  simulation  (LES)  modeling  of 
such  flows.  To  this  end,  direct  numerical  simulation  (DNS)  modeling,  relying  on  spectral 
methods  that  are  not  subject  to  dispersion  errors  and  which,  in  combination  do  not  exhibit 
unphysical  scalar  excursions,  was  also  part  of  the  effort 

2.3.1  Non-reacting  variable-density  compressible  and  incompressible  flows 

Direct  numerical  simulations  (DNSs)  that  investigate  fundamental  characteristics  of  baroclinic 
vorticity  generation  that  occurs  in  turbulent  flows  subject  to  externally  imposed  acceleration 
fields  are  in  progress.  In  the  context  of  the  AFOSR  research  focus,  such  dynamics  will  be 
encountered  when  fluid  parcels  with  density  variations  traverse  an  expansion  zone  in  supersonic 
flow,  as  seen  in  the  schlieren  image  in  Figure  1  (Section  2.1.1),  for  example,  or  are  accelerated 
through  a  contraction  in  subsonic  flow.  Very  large  (Lagrangian)  accelerations  can  be  imposed  on 
the  flow  that  are,  of  course,  indistinguishable  from  any  other  acceleration  field  as  far  as  the  flow 
dynamics  are  concerned. 

In  variable-density  flow,  interactions  between  density  and  pressure  gradients  contribute  to  the 
generation  of  baroclinic  torques  as  well  as  stabilization  (stratification)  or  destabilization 
(Rayleigh-Taylor  instability)  effects  that  have  no  counterparts  in  uniform-density  flows.  Subgrid- 
scale  (SGS)  modeling  was  developed  (Chung  and  Matheou  2014)  but  possible  turbulence 
amplification  by  unstable  configurations  remains  a  research  topic. 
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Figure  16.  Initial  density  field  and  triply-periodic  flow  domain  for  the  variable-density  flow  simulations. 


The  initial  setup  of  the  simulations  performed  is  shown  in  Figure  16.  The  evolution  of  two  fluid 
columns  with  different  densities  is  simulated  in  a  periodic  domain.  The  initially  quiescent  fluids 
move  under  the  force  of  an  external  acceleration  field,  which  is  perpendicular  to  the  density 
gradient.  This  misalignment  of  the  vertical  hydrostatic  pressure  gradient  and  horizontal  density 
gradient  induces  large-scale  baroclinic  torques  in  the  turbulent  shear-layers  that  ensue  (Figure 
17).  Small-scale  baroclinic  torques  are  also  generated  in  the  turbulent  zones  as  a  consequence 
that  would  need  to  be  accounted  for  in  a  LES. 


Figure  17.  Density  field  snapshots  on  a  vertical  (x,z) -plane  for  density  ratios  (from  left  to  right)  of 
R  =  P1/P2  =  1.4,  2,  5,  and  10,  at  time  t/r  —  0.3.  Black  marks  high-density  fluid.  Only  half  of  the 
computational  domain  is  shown.  Color  legends  in  each  panel  adjusted  to  span  the  entire  density  range. 
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Also  important  in  the  dynamics  are  instances  when  density  and  pressure  gradients  are  (nearly) 
aligned  and  parallel,  or  (nearly)  opposite.  In  such  circumstances,  the  interaction  of  the  two 
gradients  leads  to  locally  enhanced  mixing,  or  flow  stabilization,  respectively. 

Both  baroclinic -torque  generation  and  stabilization  or  destabilization  of  local  fluid  turbulence 
have  no  counterparts  in  uniform-density  flows  and  must  be  accounted  for  in  reliable  LES 
modeling. 

The  evolution  of  two  fluid  columns  with  different  densities  is  simulated  in  a  periodic  domain, 
allowing  a  full  investigation  of  baroclinic  effects.  The  flow  is  initialized  in  a  quiescent  state,  with 
free  streams  developing  relative  motion  in  an  external  acceleration  field  in  the  vertical  direction. 
Large-scale  shear  develops  with  a  turbulent  mixing  zone  between  them  (Figure  17).  Large-scale 
misalignments  of  vertical  pressure  gradients  and  horizontal  density  gradients  induce  baroclinic 
torques. 

Simulations  at  several  density  ratios  have  been  carried  out  for  1.005  <  R  =  P1/P2  —  10’  with 
numerical  grid  resolutions  up  to  10243.  A  density  ratio  of  10  poses  considerable  numerical- 
stability  challenges  and  represents  a  significant  accomplishment.  The  DNS  relies  on  a  pseudo- 
spectral  numerical  method  that  guarantees  discrete  mass  conservation,  regardless  of  iteration 
errors,  e.g.,  Chung  &  Pullin  2010).  The  numerical  calculation  imposes  a  mean  momentum  that  is 
constant  in  time,  ensuring  that  no  external  forces  other  than  gravity  drive  the  flow.  This  is  done 
by  adding  in  a  uniform  pressure  gradient,  /]  =  —  <  p  >  gSi3,  to  the  momentum  equation,  where 
<  p  >  is  the  mean  fluid  density  in  the  domain.  The  momentum  equation  then  becomes, 

DtUl=-i(aiP+ra-^3+^.  (2.3.1) 

The  term  r in  this  equation  is  the  viscous  stress  tensor,  and  Dt  =  D/D t  =  d/dt  +  (u  ■  V)  is  the 
Lagrangian  time-derivative  operator. 

The  simulations  assume  a  uniform  Schmidt  number  of  unity,  i.e.,  Sc  =  p/pD  =  1,  everywhere, 
where  D  is  the  species  diffusivity.  This  produces  comparable  viscous  and  diffusion  scales,  but, 
because  p  =  p(x,  t),  requires  a  diffusivity  that  is  a  function  of  space  and  time,  with  the  dynamic 
viscosity,  p,  assumed  uniform. 

The  simulations  verify  the  setup  and  show  a  flow  evolving  from  an  initially  homogeneous  and 
isotropic  perturbation  of  the  density  interface  (e.g.,  Cook  &  Dimotakis  2001,  Cabot  et  al.  2005) 
to  fully  turbulent  flow  with  spanwise-coherent  structures,  similar  to  uniform-density  free-shear 
layers.  Free-stream  acceleration  produces  smallest  (Kolmogorov)  scales  that  decrease  in  size  as 
time  evolves,  unlike  free-shear  layers  with  a  constant  freestream  velocity  difference.  The 
viscosity  is  chosen  for  each  DNS  along  with  its  grid  resolution  so  the  flow  remains  well 
resolved,  with  the  maximum  resolved  (de-aliased)  wavenumber  and  Kolmogorov  scale,  AK,  such 
that  kmaxAK  >  1.5,  throughout  the  flow  domain  and  during  its  evolution. 

Figure  18  shows  the  evolution  of  the  shear-/mixing-layer  width,  which  exhibits  a  laminar 
(diffusive)  evolution  up  to  a  normalized  time  t h  ~  0.2,  with  r  =  2n(l/ Jig)1/2 ,  where  l  =  L/2 
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and  L  the  spatial  extent  of  the  periodic  domain,  g  the  imposed  acceleration  magnitude,  and 
cA  =  (R-  1)/(R  +  1)  the  Atwood  number.  Following  transition  to  turbulence,  the  mixing-layer 
width  6,  defined  as  the  width  within  the  mixture-fraction  range  of  0.01  <  Y  <  0.99,  where 
(P2  <  Pi) 

1  1 

Y(x,t)  =  P\  p(X't}  (2.3.2) 

P2  Pi 

is  the  mass-fraction  of  high-density  fluid  in  the  mixture,  and  Reynolds  number  increase  with 
time  significantly  faster  following  transition  from  the  laminar-flow  regime  (Figure  18). 


Figure  18.  Non-dimensionalized  shear-layer  width,  8/1,  scaled  by  the  initial  non-dimensional 
time.  Colored  lines  are  from  simulations  and  black  lines  provide  slope  (growth-rate)  references. 

The  simulations  indicate  an  initial  diffusive  regime  with  a  shear-layer  width  growth  rate  of 
8(t)  oc  t1//2,  followed  by  a  relatively  sharp  transition  to  a  turbulent  regime,  in  which  shear-layer 
width  increases  as  the  cube  of  time,  i.e.,  8(t)  oc  t3.  Flow  slices  in  Figure  17  correspond  to  this 
turbulent-flow  regime. 

Many  turbulent  flows  grow  linearly  in  space  or  time.  Rayleigh-Taylor  (RT)  flows  grow 
quadratically  in  time,  i.e.,  dRX(t)  oc  t2  (e.g..  Cook  &  Dimotakis  2001  and  references  therein). 
This  flow  exhibits  a  new  turbulence  regime  that  is  the  fastest-growing  mixing  region 
documented  to  date  to  our  knowledge.  Similarity  and  dimensional  analysis  arguments  support 
the  unique  growth  of  the  turbulent  region  in  this  flow  environment. 

Analytical  results  verified  some  of  the  bulk  flow  statistics,  such  as  the  time-dependent  velocity 
difference  across  the  growing  shear  layers  free  stream  velocity  (Figure  19,  left),  the  mean  density 
of  the  flow  computed  using  the  theoretical  estimates  of  large-scale  structure  convection  velocity 
(Dimotakis  1986),  whose  comparison  with  the  DNS  results  is  plotted  in  Figure  19,  right),  and  the 
mixing-width  evolution  (Figure  18). 
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Figure  19.  Time  evolution  of  normalized  free  stream  vertical  velocity  (left)  and  mean  density  within  the 
mixing  region  (right).  Solid  lines  denote  numerical  results  and  dashed  lines  denote  analytical  expectation. 

Shear-layer  mixing  behavior  depends  on  the  freestream  density  ratio.  Asymmetric  entrainment  is 
characteristic  of  a  spatially  developing  shear  layer  for  uniform-density  flows.  However,  these 
simulations  are  of  temporally  developing  shear  layers  and  support  asymmetric  entrainment  as  a 
result  of  differences  in  freestream  densities  (Dimotakis  1986). 


Figure  20.  Specific  kinetic  energy  spectra,  Su.u(k)  (solid  lines),  and  kinetic  energy  spectra, 
S,.,(fc)/p  (dashed  lines),  compiled  at  the  same  Reynolds  number  of  Re  —  8500  in  all  cases, 
where  j  =  p^u  and  j  •  j  =  p|u  •  u|  =  pu2. 

A  challenge  in  variable-density  flows  has  been  how  to  resolve  kinetic  energy  into  its  spectral 
components.  For  uniform-density  flows,  the  Fourier  transform  of  the  autocorrelation  of  the  scalar 
product  of  the  velocity  field,  i.e.,  5u.u(k)  =  (F{(u(x  +  r,  t)  •  u(x,  t))x  t;  r  ^  k}  in  general, 
provides  the  wavenumber  decomposition  of  u  •  u  =  u2,  which  is  the  specific  kinetic  energy.  For 
uniform-density  flow,  this  is  simply  multiplied  by  the  uniform  density  to  yield  the  kinetic  energy. 
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The  full  kinetic  energy  in  variable-density  flow,  however,  is  given  by  (1/2)  of  pit2.  This 
challenge  was  addressed  by  considering  the  scaled  velocity  variable  j  =  p1^2  u  (Kida  and  Orzag 
1992),  whose  scalar  product,  j  •  j  =  p|u  •  u|  =  pit2  then  yields  the  desired  quantity.  This  allows 
the  spectral  decomposition  of  the  (full)  kinetic  energy  of  the  variable  density  flow  to  be 
expressed  as  Sj.j(k)  =  y{(j(x  +  r,  t)  •  j(x,  t))xt;  r  ->  k}.  This  was  then  scaled  by  the  mean 
density  of  the  fluid,  p,  i.e.,  Syi(k)/p,  so  it  can  be  directly  compared  with  the  conventional 
(specific)  kinetic  energy  spectrum.  This  comparison  is  plotted  in  Figure  20  and  attests  to  the 
success  of  this  velocity  surrogate  variable  in  bringing  variable-density  statistics  in  direct 
comparison  with  uniform-density-flow  statistics. 

These  and  other  results  have  been  presented  and  documented  in  the  paper  by  Gat  et  al.  (2016).  A 
more-detailed  paper  with  additional  results  and  discussion  is  presently  in  preparation  (Gat  et  al. 
2017). 

This  part  of  the  effort  has  been  co-sponsored  by  the  DOE  Research  Grant  Award  DE- 
NA0002382  and  continues  at  this  writing  under  its  sponsorship. 

2.3.2  Modeling  and  simulation  of  chemically  reacting  flows  and  combustion 

Cymbalist  and  Dimotakis  (2015)  and  Cymbalist  (2016)  proposed  an  approach  to  model  high¬ 
speed  combustion  flows,  including  the  effects  of  ignition  delay  and  combustion  in  distributed- 
reaction  zone  (DRZ)  regimes.  Their  evolution- variable  manifold  (EVM)  approach  explicitly 
captures  and  represents  ignition  delay  and  subsequent  heat  release  as  separate  phases  of 
combustion.  Subgrid-scale  combustion  is  represented  by  well-stirred  reactors  in  unresolved 
computational  cells  that  are  given  an  initial  condition  and  are  integrated  in  time  until  they  reach 
equilibrium.  The  time  domain  is  then  converted  to  a  (Lagrangian)  evolution  variable,  x,  that 
tracks  the  variation  of  the  gas  reactant  state  from  an  initial  condition  at  r  =  0  as  it  evolves  to 
equilibrium  at  r  =  l.1  In  the  flow  simulation,  x  is  computed  as  a  field  variable  and  its  value 
represents  the  Lagrangian  evolution  of  a  reacting  fluid  element  as  it  is  transported,  mixes,  and 
reacts  through  the  computational  domain.  The  EVM  approach  pre-computes  and  tabulates  the 
gas  state  as  a  low-order  manifold  in  terms  of  the  thermochemical  state  and  reaction  progress  in 
terms  of  energy,  density,  mixture  fraction,  and  the  evolution  variable. 

A  well-stirred  reactor  (WSR)  model  is  expected  to  represent  the  small-scale  structure  of  DRZ 
combustion  and  is  used  to  generate  the  manifold  based  on  either  direct  experimental  data  in  the 
case  of  autoignition,  or  detailed  kinetics  models,  to  the  extent  these  are  known.  With  this 
approach,  the  complete  thermochemical  state  of  the  gas  is  accessible  to  the  flow  field 
computation,  making  it  possible  to  develop  a  numerical  method  that  properly  conserves  energy 
in  both  compressible  and  incompressible  flow,  and  is  consistent  with  the  evolving 
thermochemical  state.  Because  the  low-dimensional  manifold  that  describes  the  combustion 
process  is  pre-computed,  the  EVM  approach  significantly  reduces  the  simulation  cost  relative  to 

1  The  symbol  r  in  this  section  refers  to  the  Lagrangian  time  of  a  fluid  parcel  that  is  transported  and  evolves,  and  is 

unrelated  to  the  same  symbol  in  the  previous  section  that  refers  to  the  (fixed)  time  scale  in  the  DNS  computations. 
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full  kinetics  models  and  even  relative  to  in  situ  tabulation  based  methods.  Thus,  the  EYM 
approach  has  the  potential  to  greatly  reduce  the  cost  and  improve  the  reliability  of  large  eddy 
simulation  (LES)  modeling  of  high-speed  high  Reynolds  number  turbulent-combustion  flows. 

Reliable  simulation  of  high-speed  turbulent  combustion  flows  requires  four  modeling  elements: 
compressible  flow,  near-wall  turbulent  motion,  free  or  detached  turbulent  motion,  and 
combustion  progress  including  ignition-delay  effects.  In  the  research  project,  standard 
approaches  were  used  for  the  first  three  elements,  and  focused  on  accurate  and  efficient 
modeling  of  turbulent  combustion  at  conditions  relevant  to  high-speed  flight.  To  this  end,  an 
existing  computational  fluid-dynamics  code  was  modified  and  extended  (the  University  of 
Minnesota  US3D  code)  for  wall-modeled  LES  of  compressible  reacting  flows  to  represent  the 
detailed  thermochemical  state  of  a  reacting  fuel-air  mixture,  based  on  the  reaction-evolution 
variable,  a  mixture-fraction  variable,  and  the  coupled  thermodynamic  state.  The  focus  was  on  the 
development  of  a  numerical  method  that  is  rigorously  consistent  with  the  thermochemical  model 
so  the  numerical  method  correctly  reflects  the  variations  in  the  thermodynamic  state.  We  applied 
the  proposed  numerical  flux  function  to  the  supersonic  reacting  hydrogen  injection  experiments 
of  Gamba  and  Mungal  (2015)  using  the  hydrogen-air  chemical  kinetics  model  of  Burke  et  al. 
(2012).  More  recent  simulations  have  involved  the  use  of  a  data-driven  tabulation  of  ethylene 
ignition  delay  time  along  with  tabulated  oxidation  kinetics  and  wall-modeled  large-eddy 
simulations  of  the  University  of  Virginia  dual-mode  combustor  configuration  with  ethylene  fuel. 

Developing  the  CFD  methods  so  that  the  EVM  approach  could  be  implemented  in  a  CFD  code 
involved  several  steps  and  the  use  of  a  novel  numerical  flux  function  for  high-order  of  accuracy 
simulations.  For  compressible  flows,  the  thermodynamics  and  gas  dynamics  are  tightly  coupled, 
and  this  coupling  must  be  done  in  a  rigorously  self-consistent  fashion.  With  the  EVM  approach 
(and  the  standard  flamelet-progress  variable  approach),  the  thermodynamic  state  is  obtained  from 
a  table  that  is  pre-computed  at  the  relevant  flow  conditions.  Thus,  it  was  necessary  to  modify  the 
CFD  methodology  to  use  the  tabulated  state  and  to  compute  all  necessary  thermodynamic 
variables  from  this  state. 

It  was  found  that  tabulating  the  chemical  state  as  a  function  of  the  evolution  variable  x,  the 
mixture  fraction  Z,  the  gas  density  p,  and  the  specific  energy  e,  yields  the  most  accurate  and 
general  approach.  The  rate  of  advance/production  of  x  is  also  tabulated.  With  this  approach,  all 
thermodynamic  variables  can  be  computed  from  the  tabulated  data,  along  with  the  transport 
properties  (similar  to  Oevermann  2000).  This  approach  is  much  more  accurate  and  general  than 
attempting  to  tabulate  all  required  thermodynamic  and  transport  property  variables,  as  in 
Saghafian  et  al.  (2015),  for  example. 

Secondly,  a  new  high-order  (4th-order)  convective  flux  formulation  was  developed  that  uses  the 
tabulated  information,  yet  produces  a  fully  consistent  numerical  flux  function.  Here,  the  state  of 
the  gas  at  a  cell  face  is  constructed  with  high-order  interpolants  for  species  mass  densities, 
velocity  components,  and  pressure.  Then,  the  face  value  of  temperature  is  found  from  the 
equation  of  state  of  the  local  mixed  gas  constant  and  pressure,  which  is  then  used  to  construct  the 
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face  value  for  the  mixture  energy.  Next,  the  flux  at  the  element  face  is  computed.  This  is  found 
to  be  much  more  accurate  and  robust  than  computing  the  energy  in  an  approximate  fashion. 
Details  of  this  formulation  are  documented  in  Candler,  Cymbalist,  and  Dimotakis  (2017). 

An  additional  consideration  involved  the  treatment  of  the  enthalpy  flux  attributable  to  mass 
diffusion.  It  was  found  that  this  term  must  be  computed  accurately  based  on  the  gradients  of 
individual  chemical  species,  rather  than  relying  on  approximations  involving  gradients  of  the 
evolution  variable  and  mixture  fraction.  Without  a  correct  evaluation  of  this  term,  unphysical 
bimodal  temperature  distributions  across  heat-release  fronts  may  occur. 

Extensive  testing  of  the  EYM-LES  methodology  was  conducted.  For  example,  the  use  of 
tabulated  chemical  species  data  was  mimicked  for  simple  one-step  reactions  to  verify  that  the 
EYM  code  gives  the  same  answer  as  the  original  code.  These  comparisons  were  instrumental  in 
determining  the  correct  approach  for  developing  the  numerical  methods.  As  far  as  we  know,  such 
a  method  verification  step  has  not  been  performed  for  previous  implementations  of  the  flamelet- 
progress  variable  approach,  for  example. 

Many  simulations  with  the  LES-EVM  approach  have  been  performed.  For  example,  Figure  21 
plots  the  measured  OH  planar  laser  induced  fluorescence  (OH-PLIF)  signal  on  the  center  plane 
of  the  Gamba  and  Mungal  (2015)  supersonic  jet  in  crossflow  experiments.  Also  plotted  are 
several  realizations  of  the  computed  OH-PFIF  wall-modeled  FES  using  the  EVM-FES 
method/code  developed  as  part  of  the  AFOSR-funded  research  project  described  in  this  report. 
We  have  attempted  to  match  the  color  map  used  to  plot  the  experimental  and  computational  data, 
however  the  experimental  OH-PLIF  was  not  calibrated,  so  quantitative  comparisons  are  not 
possible. 

Consistent  with  experiment,  the  EVM-LES  predicts  the  presence  of  the  thin  flame  front  at  the 
leading  edge  of  the  jet,  along  with  more  diffuse  distributed  combustion  between  the  jet  plume 
and  the  surface  downstream  of  the  injector.  Detailed  analysis  of  the  Karlovitz  number  in  this 
region  confirms  that  the  EVM  approach  is  capable  of  predicting  combustion  in  the  distributed 
reaction  (DRZ)  regime. 

Simulations  of  the  University  of  Virginia  combustor  with  ethylene  fuel  were  also  performed  and 
are  presently  being  analyzed  and  compared  to  the  available  experimental  data.  Initial  results 
appear  to  be  promising  and  we  hope  to  document  this  part  of  the  work  in  the  future. 

The  EVM-LES  method  was  demonstrated  to  be  4  times  faster  than  a  full  11 -species  chemical 
kinetics  simulation  for  hydrogen-air  combustion.  For  ethylene  combustion,  it  is  estimated  that 
EVM-LES  is  about  12  times  faster  than  using  a  skeletal  mechanism  such  as  in  the  work  of 
Potturi  and  Edwards  (2015).  Incorporating  a  detailed  full-kinetics  ethylene-air  model  that  would 
require  transport  all  participating  chemically  reacting  species  would  place  such  a  simulation 
beyond  present-day  computational  reach. 
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Figure  21.  Centerline  OH-PLIF  images  from  the  Gamba  and  Mungal  (2015)  experiment  (top)  and  3 
images  of  the  107  million  element  wall-modeled  LES-EVM  separated  by  50  ps  in  time. 

Details  of  the  integration  of  the  EVM  approach  with  the  CFD  code  for  wall-modeled  LES  are 
available  in  Candler,  Cymbalist,  and  Dimotakis  (2015,  2017),  Cymbalist  (2016),  and  Cymbalist, 
Candler,  and  Dimotakis  (2016). 

Extensive  work  was  also  conducted  on  the  large-eddy  simulation  of  two  Caltech  reacting 
hydrogen-fluorine  mixing-layer  experiments  documented  in  the  literature  (Slessor  et  al.  1998  and 
Bonanos  et  al.  2008).  Complete  successful  simulations  of  these  flows  have  not  yet  been  obtained 
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for  several  reasons.  The  work  is  continuing,  albeit  slowly,  and  we  anticipate  that  archival  papers 
will  be  submitted  on  these  simulations  in  the  near  term. 

Issues  encountered  in  the  simulation  of  these  flows  led  to  important  advancements  in  several 
areas;  here  the  focus  of  the  discussion  is  on  the  most  important  problem,  spurious  or  unphysical 
numerical  overshoots  and  undershoots  in  conserved  scalars,  as  discussed  in  detail  in  Section  2.4, 
below.  Additional  difficulties  were  encountered  with  proper  specification  of  the  outflow 
boundary  condition  so  that  the  exiting  vertical  flow  does  not  interact  with  the  flow.  This  problem 
was  addressed  by  including  a  very  large  outflow  domain  to  represent  the  experimental  outflow 
expansion  in  the  test-section  duct  and  a  portion  of  the  dump  tank,  and  the  use  of  advanced 
boundary  conditions  from  the  literature.  We  also  found  that  the  breakdown  of  the  mixing  layer  is 
extremely  sensitive  to  inflow  conditions  and  grid  resolution  in  the  vicinity  of  the  splitter  plate. 
This  is  not  surprising  since  the  initial  instabilities  that  form  near  the  inflow  drive  the  rate  of  roll¬ 
up  of  the  mixing  layer.  Unfortunately,  this  makes  the  simulations  very  difficult  because  small 
changes  to  inflow  boundary-layer  properties  can  produce  substantially  different  results.  We  have 
used  digital  filtering  approaches  to  provide  physically  meaningful  inflow  data  for  turbulent 
boundary  layers.  While  this  improved  comparisons  with  experiment,  we  are  still  studying  this 
problem  to  fully  understand  how  to  best  impose  conditions  provided  by  the  state  of  the  inflow 
boundary  layers. 

2.4  Unphysical  scalar  excursions  in  finite-difference  LES  modeling 

Unphysical  scalar  excursions  were  mentioned  in  the  discussion  above  are  responsible  for 
significant  challenges  in  the  simulation  of  turbulent  mixing,  chemically  reacting  flows,  and 
combustion.  A  discussion  of  the  issues  follows  below.  Further  details  are  documented  in  the 
publications  cited. 

2.4.1  Excursions  in  passive-scalar  transport 

The  LES  equation  for  a  resolved  passive  scalar  field,  Z(x,y,z,i),  is  given  by  (summation- 
convention  implied), 

dtZ  +  dj ( UjZ )  =  VdfjZ  -  djOj ,  (2.4.1) 

where  u (x,y,z,  t)  =  )Uj(x,y,  z,  t)  is  the  vector  flow  velocity,  V  is  the  scalar-species  diffusivity 
(here  assumed  constant  and  uniform),  dfj  is  the  Laplacian,  and  a(x,y,z,  t)  =  j<jj(x,y, z,  t)  is  the 
SGS  scalar  flux. 

Solutions  to  Eq.  (2.4.1)  must  satisfy  the  maximum  principle  i.e.,  a  Z  field  that  is  bounded  in  the 
flow  domain  by  its  initial  and  boundary  values.  In  practice,  numerical  solutions  to  Eq.  (2.4.1) 
obtained  from  difference  methods  incur  dispersion  errors  that  result  in  unphysical  scalar 
excursions,  in  violation  of  the  maximum  principle. 

Matheou  &  Dimotakis  (2014,  2016)  reported  on  the  extent  of  such  violations  incurred  with 
several  numerical  schemes  for  a  LES  modeling  of  a  temporally  evolving  turbulent  3D  shear  layer 
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at  high  Reynolds  numbers.  It  was  found  that  the  fully  conservative  schemes  of  Morinishi  et  al. 
(1998),  where  the  truncation  error  term  is  dispersive,  generate  maximum  excursions  whereas  the 
flux-limited  approach  using  the  first-order  upwind  and  the  Lax-Wendroff  method,  because  of 
their  diffusive  nature,  suppressed  excursions.  However,  diffusive  schemes  that  rely  on  numerical 
dissipation,  i.e.,  added  artificial  diffusion,  can  over-estimate  mixing,  presenting  a  serious 
challenge  to  reliable  (predictive)  numerical  simulations. 

Scalar  excursions  were  studied  in  a  triply  periodic  shear  flow  arranged  in  a  vertically  periodic 
stack  of  mixing  layers  (Figure  22).  Three  parameters  were  varied  in  the  LES  runs:  the 
discretization  of  the  convection  terms,  grid  resolution,  and  the  SGS  model  type,  including  model 
parameters.  Two  types  of  scalar  boundedness  error  were  studied:  global  and  local  (or  internal) 
unphysical  scalar  excursions.  Global  excursions  correspond  to  violations  of  the  minimum  and 
maximum  principle  with  respect  to  the  initial  and  boundary  conditions.  That  is,  any  scalar  value 
that  is  less  or  larger  than  the  boundary  values,  e.g.,  mixture  fraction  values  Z  <  0,  or  Z  >  1. 
Local  excursions  are  defined  as  violations  of  the  scalar  boundedness  property  in  the  mixed  fluid, 
i.e.,  0  <  Z  <  1. 


o 


Figure  22.  Evolution  of  the  passive  scalar  field.  Graphic  computed  for  t  =  20s.  Z  =  1  (white)  fluid 
moves  from  left  to  right  with  a  freestream  speed  of  U  =  +1  m/s.  Z  =  0  fluid  (black)  moves  from  right  to 
left  with  U  =  —  1  m/s. 

Global  excursions  are  straightforward  to  diagnose  and  compile  error  metrics  and  their 
characteristics.  They  were  analyzed  by  considering  the  minimum  and  maximum  in  the  entire 
computational  domain  and  the  volume  of  fluid  with  scalar  values  exceeding  an  excursion 
threshold.  For  a  family  of  non-dissipative  fully  conservative  schemes  (Morinishi  et  al.  1998), 
excursions  decrease  as  the  order  of  accuracy  of  the  numerical  approximation  increases.  A 
Quadratic  Upstream  Interpolation  for  Convective  Kinematics  (QUICK,  Leonard  1979),  which  is 
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a  linear  upwinding  scheme,  results  in  smaller  excursions  than  non-dissipative  schemes.  A  flux- 
limited  monotone  scheme  is  also  assessed  and,  as  expected,  no  scalar  excursions  are  observed. 
Global  unphysical  scalar-excursion  statistics  are  found  to  strongly  depend  on  the  SGS  model  and 
model  parameters. 

Two  SGS  models  were  assessed,  the  stretched-vortex  model  (SVM),  documented  in  a  series  of 
papers  by  Pullin,  Saffman  and  co-workers,  and  a  constant-coefficient  Smagorinsky  (1963) 
model.  The  excursions  are  significantly  reduced  when  the  characteristic  SGS  scale  is  set  to  twice 
the  grid  spacing  in  runs  with  the  stretched-vortex  model,  even  though  that  requires  a 
(considerably)  greater  computational  effort.  Maximum  global  unphysical  excursions  and  volume 
fraction  show  opposite  trends  with  respect  to  resolution.  The  maximum  excursion  increases  as 
resolution  increases,  whereas  the  volume  fraction  decreases,  as  might  be  expected  on  the  basis  of 
statistics;  the  number  of  grid  points  and  sample  size  increases  with  resolution.  In  contrast,  the 
volume  fraction  of  unphysical  excursions  decreases  with  resolution  because  the  SGS  models 
perform  better  at  higher  grid  resolutions. 

Local-excursion  characteristics  are  defined  with  respect  to  the  local  non-diffusive  limit.  To 
define  a  local  excursion,  the  scalar-transport  equation  is  integrated  in  time  from  an  initial  state,  t, 
for  a  short  time  interval  At,  i.e., 

rt+At  rt+At  rt+At 

J  dtZdt  +  J  di(uiZ)dt  =  Vj  d^Z  dt.  (2.4.2) 

The  value  of  Z  at  t  +  At  depends  on  contributions  to  Zt  from  convection  and  diffusion,  i.e., 

rt+At  rt+At 

Zt+At  =  zt-  J  di(utZ)  dt+Vj  d\Z  dt.  (2.4.3) 


In  the  limit  of  vanishing  positive  diffusivity,  D  0+,  the  sign  of  the  diffusive  term  determines 
the  relation  of  Zt+At  with  respect  to  the  non-diffusive  limit.  Thus,  the  local  constraint  of  the 
Eulerian  Z  is. 


't+At  >  sign 


(rt+At 
1  ‘ 


df,Zdt 


t+At 


dj.  (ttjZ)dt  —  -^at.nd  > 


(2.4.4) 


where  the  right-hand  side  denotes  the  non-diffusive  short-time  estimate.  This  constraint  is  local 
in  space  and  time  and  results  in  the  requirement  that  only  local  mixing/diffusion  of  Z  is  allowed. 
That  is,  it  requires  that  local  Z  values  must  be  on  the  “convex  side”  of  the  non-diffusive  limit. 

To  diagnose  (2.4.4)  in  LES,  the  convection  term  must  be  computed  sufficiently  accurately  for  an 
interval  At,  which  requires  performing  an  inviscid  integration.  A  more-accurate  inviscid 
integration  for  the  scalar  is  performed  on  a  finer  grid.  In  the  present  runs,  we  utilize  a  refinement 
factor  of  4,  i.e.,  Axfine  =  Ax/4,  for  diagnostic  purposes. 
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At  time  t,  velocity  and  scalar  fields  are  interpolated  on  the  fine  grid  using  trilinear  interpolation, 
which  preserves  total  mass.  The  scalar  equation  on  the  fine  grid  is  then  integrated  for  time  At 
using  the  current  fine-grid-interpolated  velocity  field  at  each  sub-step  of  the  Runge-Kutta.  This 
process  results  in  the  RHS  of  (2.4.4),  i.e.,  ZAtnd.  Constraint  (2.4.4)  is  then  evaluated  by 
comparing  the  co-located  ZAtnd  on  the  fine  grid  and  Z  on  the  original  (coarse)  LES  grid.  No 
averaging  is  performed  from  fine  to  coarse  grids,  because  this  would  result  in  diffusion  of  the 
fine-grid  Z-field  and  invalidate  the  inviscid  nature  of  the  integration  that  is  required. 

Figure  23  shows  the  frequency  distribution  (i.e.,  the  two-dimensional  histogram)  of  the  local 
excursion-scalar  gradient  pairs  for  11  flow  realizations  in  the  t  =  10-20s  interval,  for  a  run 
with  2563  grid  points,  fourth-order  convection,  and  the  stretched-vortex  model.  Unphysical 
excursions  (y-axis)  are  relatively  rare  and  about  97%  of  the  occurrences  are  accounted  for  in  the 
y  =  0  bins  in  the  histograms  of  Figure  23,  corresponding  to  the  dark  band  along  the  x-axis.  The 
pairing/association  of  excursions  to  scalar  gradient  is  done  in  two  ways:  using  co-located  pairs 
(Figure  23,  left  panel)  and  by  pairing  each  excursion  value  with  the  maximum  gradient 
magnitude  in  the  surrounding  7x7x7  grid  cells  (Figure  23,  right  panel).  This  is  to  account  for  the 
offset  of  the  maximum  dispersion  error  with  respect  to  the  local  gradient,  which  is  usually  2-3 
grid  cells.  The  joint  histograms  in  Figure  23  show  that  excursions  and  scalar  gradient  magnitude 
are  not  well  correlated.  Even  when  accounting  for  non-colocation,  maximum  excursions  tend  to 
occur  at  intermediate  values  of  scalar-gradient  magnitude.  Moreover,  for  a  given  scalar  excursion 
value,  the  gradient  distribution  is  broad. 


Figure  23.  Two-dimensional  frequency  distributions  of  local  excursion-scalar  gradient  pairs  for  11  flow 
realizations  in  t  —  10- 20s  for  a  run  with  2563  grid  points,  fourth  order-convection  and  the  stretched- 
vortex  model.  Darker  shades  correspond  to  higher  frequency.  Frequency  values  span  a  large  range  of 
scales  and  contours  are  logarithmically  spaced,  i.e.,  log(/i  +  1),  where  h  is  the  frequency  of  occurrence 
of  a  given  pair.  Left  panel  shows  co-located  pairs.  Right  panel  shows  excursions  paired  with  the 
maximum  scalar  gradient  magnitude  in  surrounding  7x7x7  grid  cells. 
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Figure  24  shows  the  frequency  distribution  of  the  local  excursions  and  estimate  of  the  scalar- 
dissipation  rate.  Unlike  the  local-excursion  and  scalar-gradient  pairs  in  Figure  23,  a  significant 
fraction  of  local  unphysical  excursions  can  be  seen  to  occur  near  grid  points  with  low  scalar 
dissipation  values.  There  is  a  very  small  fraction  of  grid  points  with  (unphysical)  negative 
dissipation  rates  but  for  most  of  these  points  no  unphysical  excursions  are  observed  (contour 
scales  in  Figure  23  and  Figure  24  are  with  a  logarithmic  spacing).  The  right  panel  of  Figure  24 
shows  pairs  of  scalar  excursions  with  the  minimum  dissipation  in  the  surrounding  7x7x7  grid 
cells  (similar  to  the  right  panel  of  Figure  23).  All  grid  cells  with  unphysical  scalar  excursions  are 
at  a  distance  less  than  4  grid  cells  from  locations  with  near-zero  scalar-dissipation  rate. 

The  analysis  of  scalar-excursion  statistics  shows  that  unphysical  scalar  excursions  in  LES  result 
from  dispersive  oscillations  of  the  convection-term  discretization  at  times  and  locations  where 
the  subgrid-scale  model  provides  insufficient  dissipation  to  produce  a  sufficiently  smooth  scalar 
field.  Further  details  are  provided  in  Matheou  &  Dimotakis  (2016).  A  conclusion  of  that  work  is 
that  further  improvements  in  SGS  models  can  mitigate  and  may  be  required  to  address 
unphysical  scalar  excursions. 
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Figure  24.  Two-dimensional  distributions  of  local  excursion  and  scalar-dissipation-rate  pairs  for  11 
flow  realizations  in  t  =  10- 20s  for  runs  with  2563  grid  points,  4th-order  convection  and  stretched- 
vortex  model.  Darker  shades  correspond  to  higher  frequencies.  Contour  levels  are  the  same  as  in  Figure 
23.  Frequency  values  span  a  large  range  of  scales  and  contours  are  logarithmically  spaced,  i.e.,  log(/i  + 
1),  where  h  is  the  frequency  of  occurrence  of  a  given  pair.  Left  panel  shows  co-located  pairs.  Right  panel 
shows  excursions  paired  with  the  minimum  scalar-dissipation  rate  in  the  surrounding  7x7x7  grid  cells. 

The  causes  leading  to  insufficient  SGS  scalar  dissipation  are  currently  not  well  understood. 
Moreover,  it  is  not  clear  if  the  problem  originates  from  the  momentum  SGS  modeling  or  only 
affects  the  scalar  fields,  implying  a  breakdown  in  assumptions  that  relate  the  turbulent 
momentum  and  scalar  fields.  Even  though  the  scalar  excursion  volumetric  error  is  small,  its 
nature  can  lead  to  compounded  errors  in  LES  involving  active  scalars.  Higher  order  non- 
dissipative  finite-difference  discretizations  do  not  yield  fast  error  reduction,  while  low-order 
discretizations  produce  excessively  high  errors,  and  general  dissipative  schemes  must  be 
regarded  as  a  blunt  tool,  since  they  introduce  artificial  dissipation  in  locations  where  the  SGS 
model  may  already  provide  the  necessary  dissipation. 
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2.4.2  Unphysical  excursions  in  active-scalar  transport 

In  an  extension  of  this  research  to  active  scalars,  it  was  noted  that  unphysical  scalar  excursions 
were  responsible  for  the  major  bottleneck  associated  with  the  simulation  of  the  Caltech  HF 
experiments  (Subbareddy  et  al.  2014;  Ferrero  et  al.  2013;  Kartha  et  al.  2014,  2015).  This  has 
relevance  to  many  reacting-flow  simulations  where  unphysical  undershoots  and  overshoots  may 
occur  in  species  mass  fractions  or  other  scalars  that  must  be  rigorously  bounded.  This  is  a 
particular  problem  for  compressible  LES  where  interfaces  between  different  gas  types  are 
usually  poorly  resolved  by  the  grid.  The  high-order  (fourth-order,  or  higher)  methods  required 
for  LES  often  produce  spurious  values  of  scalars  across  these  interfaces.  During  the  research 
project,  we  worked  to  eliminate  this  type  of  numerical  error;  we  have  now  developed  and  tested 
a  new  limiting  approach  that  allows  spatially  fourth-  and  sixth-order  accurate  solutions  to  be 
obtained.  Two  papers  are  nearing  completion  on  this  work. 

Physical  bounds  of  conserved  scalar  values  set  by  the  initial  and  boundary  values  of  a  simulation 
are  not  violated.  For  instance,  the  numerical  method  must  guarantee  that  the  mass  fractions  of 
chemical  species  are  individually  in  the  range  [0,1]  and  sum  to  unity.  In  a  reacting  gas,  species 
mass  fractions  must  also  be  bounded,  but  their  limiting  values  are  not  set  by  the  boundary 
conditions.  This  makes  it  much  more  difficult  to  control  spurious  mass-fraction  values.  Naively 
correcting  deviations  generally  results  in  a  scheme  that  does  not  conserve  mass. 

The  use  of  spatially  high-order  methods  in  reacting  flow  calculations  often  results  in  local 
undershoots  and  overshoots  in  the  species  mass  fractions,  causing  unphysical  results.  This  results 
in  mass-fraction  or  other  scalar  non-conservation.  For  reacting  cases,  these  errors  can  cause 
temperature  to  exceed  the  adiabatic-flame  temperature  or  violate  other  physical  constraints. 

We  have  developed  a  switched,  low-dissipation  numerical-flux  methodology  suitable  for  scalar 
conservation  laws  in  compressible  flow.  For  passive  scalars,  the  approach  uses  a  decoupled  flux 
that  combines  several  ideas  from  the  literature  and  produces  bounded  results  with  acceptable 
levels  of  numerical  diffusion.  For  the  species  fluxes  (active  scalars),  we  have  developed  an 
efficient  novel  low-dissipation  method.  A  few  heuristic  choices,  involving  the  low-dissipation 
switch  and  the  jump  detectors  are  made:  these  parameters  can  be  tuned  for  specific  applications. 

Figure  25  presents  a  comparison  for  two  sixth-order  accurate  flux  methods  for  a  one-dimensional 
test  case  involving  large  variations  in  the  thermodynamic  properties  and  the  mass  fractions.  Note 
that  the  improved  numerical  flux  function  eliminates  the  spurious  variations  in  the  plotted 
species  mass  fraction  that  pervades  and  corrupts  the  solution  for  the  original  numerical  flux. 

Additional  multi-dimensional  test  cases  show  that  the  proposed  approach  is  also  valid  for  more 
complex  flows.  To  evaluate  the  performance  of  the  new  method  for  cases  with  variations  in 
mass-fraction  and  thermodynamic  properties,  we  perform  test  simulations  reported  in  literature. 
These  cases  include  a  two-dimensional  temporal  mixing  layer  with  multiple  species  and  a 
shock/density-bubble  interaction  that  shows  very  good  comparison  with  literature  results.  Thus, 
we  have  found  that  the  method  produces  results  with  well-behaved  variations  in  the  species- 
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densities  and  temperature,  and  the  results  improve  with  increasing  grid  refinement. 


Figure  25.  Comparison  between  the  original  and  improved  6th-order  accurate  fluxes  for  a  one¬ 
dimensional  flow  with  large  density,  temperature  and  mass  fraction  variations. 

To  further  evaluate  the  new  scalar-limiting  method,  we  performed  three-dimensional  simulations 
of  temporal  mixing  layers  with  multiple  species  (as  mentioned  above).  Here,  the  overshoots  in 
species  mass-fractions  are  mitigated  without  significantly  affecting  the  growth  of  the  mixing 
layer  and  Reynolds  stresses.  Finally,  we  investigated  high-Reynolds  number  chemically  reacting 
HF  flows  by  performing  large  eddy  simulations  corresponding  to  the  experiments  of  Slessor  et 
al.  (1998)  at  Caltech.  Example  results  from  this  simulation  are  shown  in  Figure  26.  The  Q- 
criterion  plot  shows  rich  eddy  content  indicative  of  low  levels  of  numerical  dissipation. 

Temperature  profiles  from  the  simulation,  although  close  to  the  experimental  results,  do  not 
show  satisfactory  agreement  at  present.  Again,  the  computed  temperature  distribution  is  very 
sensitive  to  the  specification  of  the  inflow  conditions.  Further  investigations  are  being  conducted 
to  address  this  issue,  and  will  be  the  focus  of  research  in  the  future. 

We  conclude  this  section  by  mentioning  a  topic  that  was  explored  but  which  is  not  part  of  the 
focus  of  the  AFOSR-funded  research,  but  one  that  benefitted  from  it  and  provides  a  rich  testing 
and  validation  environment  for  these  methods  and  considerations.  In  particular,  an  exploration  of 
active  scalar  transport  was  performed  by  exploiting  techniques  and  progress  made  as  part  of 
Sections  2.4.1  and  2.4.2  of  this  report,  focusing  on  the  consequences  of  subgrid-scale  modeling 
and  active  scalar  bound  control  in  atmospheric  dynamics.  That  work  relies  on  progress  and 
advances  documented  recently  on  flows  where  SGS  stability  corrections  have  proven  necessary 
(Chung  and  Matheou  2014),  was  co-sponsored  by  the  DOE  grant  mentioned  earlier,  published  by 
Matheou,  Chung,  and  Teixeira  (2016),  and  listed  in  the  reference  list  of  this  report  but  not  as  a 
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publication  sponsored  by  the  AFOSR  grant.  It  demonstrates  the  importance  of  observing  scalar 
bounds,  especially  active-scalar  bounds,  as  well  as  the  benefits  of  improved  SGS  modeling. 


Y 


Figure  26.  Three-dimensional  reacting  flow  simulation  a)  Q-criterion  plot  and  b)  comparison  of  present 
temperature  profile  with  experimental  data. 

2.5  Compressible  and  variable-density  SGS  modeling 

Compressible  and  variable-density  turbulence  is  described  by  non-linear  partial  differential 
equations  that  represent  the  conservation  of  mass,  momentum,  and  energy.  To  derive  the 
governing  equations  for  large-eddy  simulations,  these  equations  are  low-pass  filtered  on  a  space- 
time  volume.  Typically,  this  volume  is  related  to  the  local  grid  element  and  has  a  characteristic 
dimension.  The  filtered  variables  become  the  solution  variables,  and  there  are  unclosed  terms  in 
the  governing  equations  due  to  filtering  non-linear  quantities.  The  effects  of  these  subgrid-scale 
(SGS)  terms  cannot  be  resolved  by  the  grid,  and  must  be  represented  with  a  model. 

Almost  all  CFD  codes  for  variable-density  flows  rely  on  density-weighted  (or  Favre)  filtering. 
The  solution  variables  are,  for  example,  tq  =  piq/p,  where  the  overbar  represents  the  filtering 
operation,  the  tilde  represents  the  Favre-filtered  quantity,  p  is  the  density,  and  rq  is  the  velocity 
in  the  i  direction.  This  simplifies  the  derivation  of  the  filtered  governing  equations.  However,  it 
obscures  the  interaction  between  density  variations  and  turbulent  motion.  For  example,  the 
filtered  momentum  flux  in  the  i  direction  attributable  to  velocity  in  the  j  direction  is  piqiq,  and 
its  unresolved  effects  are  modeled  as  Ttj  =  ptqtq  —  ptqtq,  where  p,  tq,  and  tq  are  the  resolved 
solution  variables,  and  r is  (here)  the  model  for  the  unresolved  subgrid-scale  momentum  flux. 

The  model  must  account  for  the  difference  between  the  filtered  triple  product  of  the  density  and 
velocity  fields  and  the  resolved  momentum  flux.  This  is  a  complicated  quantity  that  embodies 
critical  variable-density  flow  physics.  However,  its  mathematical  form  provides  no  guidance  as 
to  how  it  should  be  represented  in  terms  of  the  resolved-scale  flow  quantities. 
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Current  LES  models  ignore  the  role  of  density  variations  in  this  product,  and  it  is  known  that  this 
approach  cannot  and  does  not  capture  the  full  dynamics  of  variable-density  turbulence. 
Furthermore,  filtering  the  energy  equation  results  in  a  quadruple  product  for  the  unresolved 
kinetic  energy  flux;  this  term  is  also  not  modeled  in  existing  CFD  codes. 

Under  the  present  project,  we  derived  an  alternate  mathematical  form  of  the  filtered  governing 
equations  for  LES  of  variable-density  turbulence  (Sidharth,  Kartha,  and  Candler  2016).  These 
developments  were  based  on  prior  studies  of  the  effects  of  SGS  baroclinic  torque  in  the  context 
of  a  variable-density  extension  to  the  SVM  (GS,  Candler,  and  Dimotakis  2014;  GS  and  Candler 
2015).  The  key  idea  is  to  apply  unweighted  filtering  (or  Reynolds-filtering),  and  then  to  use  the 
transport  equation  for  the  subgrid-scale  mass  flux,  pul  —  put,  to  eliminate  terms.  This  results  in  a 
new  form  of  the  governing  equations  with  several  advantages: 

•  The  governing  equations  are  mathematically  exact  and  represent  the  identical  physics  as 
the  conventional  Favre-filtered  equations. 

•  The  unresolved  subgrid-scale  terms  involve  only  filtered  double  products,  greatly 
clarifying  the  connection  between  the  mathematical  description  and  the  turbulence. 

•  These  terms  can  represent  subgrid-scale  counter-gradient  diffusion  and  backscatter 
(small-scale  motion  driving  motion  at  larger  scales). 

•  There  are  terms  for  the  subgrid-scale  baroclinic  torque  (due  to  misalignment  of  density 
and  pressure  gradients)  and  acceleration  (not  represented  in  current  LES  models). 

•  The  formulation  is  readily  extensible  to  reacting  flows  and  other  complex  variable- 
density  turbulent  flows. 


The  subgrid-scale  terms  are:  the  mass  flux,  specific  stress,  dilatational  flux,  pressure  work, 
pressure  acceleration,  and  internal  energy  flux.  These  terms  are  embedded  in  the  Favre-filtered 
equations,  but  it  is  impossible  to  isolate  them  in  the  filtered  triple  and  quadruple  products. 


Figure  27.  Filtered-velocity  LES  of  a  recirculating  mixing  layer:  isosurface  of  Q  criterion,  colored  by  the 
chemical-product  mass  fraction.  A  gradient  model  was  used  for  SGS  pressure-acceleration  modeling. 

We  have  carried  out  direct  numerical  simulations  of  decaying  variable-density  turbulent  flow  to 
understand  the  effect  of  variable-density  subgrid-scale  terms.  A  budget  analysis  indicates  that  the 
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new  Reynolds-filtered  subgrid-scale  terms  are  non-negligible  in  the  mean  sense  and  their 
magnitude  can  locally  exceed  the  magnitude  of  the  conventional  uniform-density  stress  term. 
The  DNS  cases  are  used  to  extract  coefficients  for  a  gradient-based  SGS  model.  This  new  SGS 
model  was  applied  to  complex  flows.  For  example,  a  recirculating  mixing  layer  is  simulated  with 
the  stretched-vortex  model  and  this  preliminary  set  of  gradient  model  coefficients  (the  flow 
configuration  and  conditions  correspond  to  the  work  of  Bonanos  et  al.  2008).  The  LES  results 
(Figure  27)  are  promising  and  are  currently  being  further  validated.  A  paper  reporting  on  this 
approach  is  in  preparation  to  be  submitted  for  publication. 
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3.  Summary  and  conclusions 

This  report  outlines  and  summarizes  research  progress  and  advances  as  a  result  of  work  under 
sponsorship  of  this  Grant,  as  well  as  the  co-sponsorship  of  some  of  the  research  topics  by  other 
funding,  as  outlined  in  Section  1 . 

In  experimental  work,  the  scaling  of  the  delay  length  for  initiation  of  chemical  reactions 
downstream  of  the  jet-injection  station  based  on  jet  Reynolds  number  was  established.  This  was 
confirmed  to  be  dominated  by  turbulent  mixing  and  not  chemical-kinetic  autoignition  delay  for 
the  flows  investigated.  A  30°  jet  inclination  dramatically  decreases  jet-injection  disturbances  to 
the  freestream,  relative  to  normal-jet  injection,  but  increases  the  distance  to  ignite,  as  a  result  of 
lower  mixing  rates.  However,  this  penalty  may  be  more  than  acceptable  considering  the  (much) 
lower  total-pressure  losses  as  evidenced  by  the  absence  of  (near-)normal  shocks  that  normal 
injection  generates.  New  experimental  techniques  were  developed  that  allowed  the  non-intrusive 
measurement  of  convecting  velocity  fields,  based  on  schlieren-imaging  velocimetry  (SICV).  The 
techniques  was  used  for  a  first  validation  of  LES  of  jet-injection  flow  in  a  supersonic  stream 
exploiting  new  geometric-optics  modeling  also  developed  as  part  of  this  grant.  Small-scale 
experiments  revealed  the  3D  and  4D  behavior  of  scalar  mixing  providing  guidance  to  subgrid- 
scale  (SGS)  modeling. 

In  large-scale  direct-numerical  simulations  (DNS)  turbulent  behavior  of  the  mixing  of  fluids  with 
small  to  high  density  ratios  subject  to  imposed  acceleration  fields  revealed  important  new  scaling 
of  various  statistics,  also  providing  guidance  for  future  subgrid-scale  modeling  of  variable- 
density  flows. 

A  new  subgrid-scale  combustion  model  was  developed  that  also  provides  for  efficient 
computation  of  complex  chemical  kinetics  in  a  turbulent-mixing  environment.  Dubbed  the 
evolution-variable  manifold  (EVM),  it  was  integrated  in  a  LES  of  turbulent  combustion  of  a 
hydrogen  jet  in  supersonic  air  crossflow,  capturing  many  important  flow  features,  perhaps  for  the 
first  time.  This  modeling  is  best  suited  for  capturing  hydrocarbon  combustion  in  the  distributed- 
reaction  zone  (DRZ)  regime  but  can  also  capture  flame  sheets  with  no  changes  in  its  structure. 

The  important  topic  of  unphysical  scalar  excursions  observed  in  finite-difference  LES  modeling 
received  considerable  attention.  It  was  possible  to  attribute  this  behavior  to  dispersion  errors  of 
finite-difference  operators.  First  attempts  were  made  to  mitigate  these  errors  that  must  be 
controlled  if  passive  and  active  scalar  transport  is  to  be  simulated  reliably  and,  in  particular,  if 
chemically  reacting  flows  are  to  be  reliably  computed. 

First  attempts  at  variable-density  and  compressible-flow  subgrid-scale  model  developments  were 
made  and  preliminary  simulations  to  assess  their  efficacy  were  made.  However,  this  topic 
remains  an  open  research  challenge  at  this  time. 

A  large  number  of  presentations,  conference  papers,  and  publications  derived  from  this  effort,  as 
listed  in  the  Reference  list  in  Section  5  of  this  report. 
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